   function [vdcom,vcvth,sdevf,ACzero,ACone,vcshocks,auone,vdech,vdech_st,vdcom_st,std_st]=vardecom_mode(G1,HH,SDX,zmat,nforc)
%  function [vdcom,vcfull,sdevf,AConevcfull,vcshocks,auone,vdech,vdech_st]=vardecom_mode(G1,HH,SDX,zmat,nforc);
%
% Compute the VCV, First order autocorrelations, Variance Decomposition and
% use with mode_evaluation.m 
%
% Inputs 
% ------
% G1    from gensys.m 
% HH    impact from gensys.m  output 
% SDX   Variance-covariance matrix implied by the parameter
%       ( from model file )
% zmat  Matrix mapping model variables to observables 
%       ( from model file )
% 
% Output  
% -------
% Full decomposition versus diagonal decompositions 
% =================================================
% Code works by computing the VCV ( nz x nz ) for each shock 
% results is an nz x nz x nx matrix  vshocks 
% Adding across pages produces nz x nz matrix  vcfull...
%
% Then working with the diagonal elements of vchocks ONLY, it is possible to
% make a variance decomposition 
% 
% Same is true of the AC 
% auone computes the full first order covariance matrix which is nz x nz x
% nx 
% then this is combined with the VCV to obtain the aggregate CORR(0) and
% CORR(1) 
%
%========================================================================== 
% vdcom (NZ x NX ) 
%            Each row  corresponds to a series, 
%            Each column to the cointribution per shocks 
%
% ACzero (NZ x NZ) Correlation Matrix for the Observable Variables 
%
% ACone  (NZ x NZ) First order Correlation matrix for Observable Variables 
%
% vcfull    Full VCV matrix 
%
% vcshocks  Variance covariance matrix, each page corresponds to a
%           shock vcfull=sum( vcshocks, 3) 
%
% auone     First order covariance matrix, ach page corresponds to a
%           shock
%
% VDECH     Variance decomposition at different horizons for observables 
%           Pass input NFORC 
% VDECH_st  Same as VDECH but for all states in the model 
% =========================================================================
%
%
% Output from gensys is of the form 
% y(t) = G1*y(t-1) + HH*e(t) 
% where v( e(t) ) = SDX (diagonal usually) 
% let  V(0) = var ( y (t) ) 
% then it follows that the autocovariance of order j is 
% V(j) = (G1^j)*V(0) 
% Also compute the corrlation and first order 
% autocorrelation matrix 
%
% Intermediate between vdecom and vdecom_simp.m 
% Need to Consolidate codes! 
% AJ 5/25/05 
% Modified June 27 2007 Work with LYAP Solvers, also added VDECH_ST 
% That reports the horizon decomposition for ALL states 
% =======================================================================
NY=size(G1,1);NX=size(SDX,1);NZ=size(zmat,1); 
%pshat=lyapunov_symm(G1,HH*(SDX'*SDX)*HH'); 
pshat=disclyap_fast(G1,HH*(SDX'*SDX)*HH'); 
% if max(max(abs(pshatc-pshat))) > 1e-2 
%     error('Numerical Discrepancies LYAPUNOV SYMM and DISCLYAP_FAST') 
% end 
pshat=0.5*(pshat+pshat'); 
deno_st=diag(pshat); 
vcvth=zmat*pshat*zmat'; 
ACone=zmat*G1*pshat*zmat';  % First order autocovariance
varf=diag(vcvth);  
sdevf=sqrt(varf); 
% Autocorrelation Correlation matrix 
% -------------------------------------
ACzero = vcvth ./ ( sdevf*sdevf'   ); 
%First order autocorrelation 
%----------------------------
ACone = ACone ./ (sdevf*sdevf' );


% Matric full state 
std_st=sqrt( diag(pshat) ); 
ac0_st=pshat ./ (std_st*std_st') ; 


vcshocks=zeros(NZ,NZ,NX); 
auone=zeros(NZ,NZ,NX); 
vdcom=zeros(NZ,NX); 
vdcom_st=zeros(NY,NX); 
for ii=1:NX; 
    pshat=real(disclyap_fast(G1,HH*((SDX(ii,:)' )*SDX(ii,:))*HH'));  
    pshat=0.5*(pshat+pshat'); 
    temp_var=zmat*pshat*zmat'; 
    vdcom_st(:,ii)=(diag(pshat))./(deno_st); 
    vcshocks(:,:,ii)= temp_var; 
    auone(:,:,ii) = zmat*G1*pshat*zmat';  % First order autocovariance     
    tempvd = diag( temp_var )./ varf ; 
    vdcom(:,ii) =tempvd  ;     
end;     
% ================================================================
%        New: Compute Forecast errors at different Horizons      %
% ================================================================
if nargin < 5 || nforc==0 || isempty(nforc);
    vdech=[];vdech_st=[];return; 
end; 
% =============================================
% 1 step ahead 
% ===============================================
Gnew=eye(NY); 
vdech=zeros(NZ,NX,nforc); 
vdech_st=zeros(NY,NX,nforc); 

vzero=HH*(SDX')*SDX*(HH'); 
vcum=zeros(NZ,NZ,nforc,NX);
vcum_st=zeros(NY,NY,nforc,NX); 
for ii=1:nforc;
    % ====================
    % Whole Variance
    % ====================
    if ii==1;
        vcfull_st=vzero;
        vcfull=zmat*vcfull_st*(zmat');
    else
        Gnew=G1*Gnew ;
        vcfull_st=Gnew*vzero*(Gnew')+vcfull_st;
        vcfull=zmat*Gnew*vzero*(Gnew')*(zmat')+vcfull;
    end;
    deno=diag(vcfull);
    deno_st=diag(vcfull_st);
    if ii==1;
        for jj=1:NX
            et=( SDX(jj,:)' )*SDX(jj,:);
            vcum_st(:,:,1,jj)=HH*et*(HH');
            vcum(:,:,1,jj)=zmat*(vcum_st(:,:,1,jj))*(zmat');
            vdech(:,jj,ii)=(diag(vcum(:,:,1,jj)))./deno;
            vdech_st(:,jj,ii)=(diag(vcum_st(:,:,1,jj)))./deno_st;
        end;
    else
        for jj=1:NX;
            et=( SDX(jj,:)' )*SDX(jj,:);
            vcum(:,:,ii,jj)=zmat*(Gnew*HH*et*(HH')*(Gnew'))*(zmat')+vcum(:,:,ii-1,jj);
            vcum_st(:,:,ii,jj)=(Gnew*HH*et*(HH')*(Gnew'))+vcum_st(:,:,ii-1,jj);
            vdech(:,jj,ii)=(diag(vcum(:,:,ii,jj)))./deno;
            vdech_st(:,jj,ii)=(diag(vcum_st(:,:,ii,jj)))./deno_st;
        end
    end
    if any( sum( squeeze(vdech(:,:,ii)) , 2 ) < 0.99 )
        error('Variance Decompositions Observables are innacurate');
    end
    if any( sum( squeeze(vdech_st(:,:,ii)) , 2 ) < 0.99 )
        error('Variance Decompositions States are innacurate');
    end
end